Transforming Gaussian diffusion into fractional, a generalized law of large numbers 
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The fractional Fokker-Planck equation (FFPE) [R. Met- 
zler, E. Barkai, J. Klafter Phys. Rev. Lett. 82, 3563 (1999)] 
describes an anomalous sub diffusive behavior of a particle 
in an external force field. In this paper we present the solu- 
tion of the FFPE in terms of an integral transformation. The 
transformation maps the solution of ordinary Fokker-Planck 
equation onto the solution of the FFPE. We investigate in 
detail the force free particle and the particle in uniform and 
harmonic fields. The meaning of the transformation is ex- 
plained based on the asymptotic solution of the continuous 
time random walk (CTRW). We also find an exact solution of 
the CTRW and compare the CTRW result with the integral 
solution of the FFPE for the force free case. 
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I. INTRODUCTION 

The continuous time random walk (CTRW), intro- 
duced by MontroU and Weiss describes different 
types of diffusion processes including the standard Gaus- 
sian diffusion, sub diffusion and Levy walks [|6|-p^. 
The CTRW has been a useful tool in many diverse fields 
and over the last three decades e.g., it was used to de- 
scribe transport in disorder medium J0,|l^,Q and in low 
dimensional chaotic systems [p|Jl5|-[l7| . More recently 
fractional kinetic equations were investigated as a tool de- 
scribing phenomenologically anomalous diffusion p8|-p4[ . 
In these equations fractional derivatives ||2^ replace or- 
dinary integer derivatives in the standard integer kinetic 
equation. It is believed that fractional calculus can be 
used to model non integer diffusion phenomena. Schnei- 
der and Wyss have formulated the fractional diffusion 
equation (see Sec. Ill) and under certain condi- 
tions |p7|-p9t the fractional diffusion equation describes 
the asymptotic large time behavior of the decoupled sub- 



diffusive CTRW (with (r^) - t^; 5 < 1). Thus some 
what like ordinary random walks which can be described 
asymptoticly by the diffusion equation, so can the CTRW 
be described by fractional diffusion equation. 

Recently, Metzler, Barkai and Klafter have inves- 
tigated fractional Fokker-Planck equation (FFPE) to de- 
scribe an anomalous sub-diffusion motion in an external 
non-linear force field. In the absence of the external force 
field the FFPE reduces to the fractional diffusion equa- 



tion. The FFPE is an asymptotic equation which extends 
the CTRW to include the effect of an external force field. 
The CTRW itself was not designed for a description of a 
particle in external force field thus the FFPE describes 
new type of behavior not explored in depth yet. 

The main result in this paper concerns a transforma- 
tion of ordinary Gaussian diffusion into fractional diffu- 
sion |^,^l|-Q. Let Pa{Y,t) be the propagator (Green's 
function) of fractional Fokker-Planck equation (a special 
case is the fractional diffusion equation), here < a < 1 
is the fractional exponent and a = 1 is the standard case 
described by ordinary Fokker-Planck equation. Then the 
solution, ^^(r,^), is found based on the transformation 





'dN{s,ty 


'a 


ds 



Pi(r,s)ds, 



and 



N{.s,t) = 1 - La 



51/C 



(1) 



(2) 



denotes the "inverse" one sided Levy stable distribution 
(i.e., La{x) is the one sided Levy stable distribution). 
For example consider the force free case. The solution of 
the FFPE Pct(r, i), with initial conditions concentrated 
on the origin, is found by transforming Pi(r, s) and the 
transformed function is the well known Gaussian solu- 
tion of the integer diffusion equation with initial con- 
ditions concentrated on the origin. The transformation 
Eq. (|l|) besides its practical value for finding solutions of 
the FFPE also explains its meaning and relation to the 
CTRW (see details below). 

We list previous work on the integral transformation 
Eq. (I) 

1 Bouchaud and Georges Q for the asymptotic long time 
limit of the CTRW, then Pi(r, s) is Gaussian, 

2 Zumofen and Klafter ||l| for CTRW on a lattice, then 
Pi(r, s) is solution of an ordinary random walk on a lat- 
tice, 



3 Saichev and Zaslavsky 1 32 1 for one dimensional solution 
of fractional diffusion equation, they have also considered 
an extension which includes Levy flights, 

4 Barkai and Silbey |Q for the fractional Ornstein- 
Uhlenbeck process. 

We generalize these results for the dynamics described 
by the FFPE, and investigate in greater detail Eq. (|l|) in 
the context of CTRW and fractional diffusion equation 
in dimensions d= 1,2,3. 
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This paper is organized as follows. In Sec. (||) we 
follow ||] and derive the transformation Eq. (P based 
on the CTRW. We find an exact solution to the CTRW 
valid for short and long times and compare this solution 
with the result obtained from Eq. (|l|). In Sec. (Ill) we 
consider the fractional diffusion equation. We investigate 
its integral solution Eq. (^ as well as the Fox function 
solution found in j2^. In Sec. (IV) we show how to 
obtain solutions of the FFPE using Eq. (0). We give 
as examples the motion of a particle in a uniform and 
harmonic force fields. In Sec. (^) a brief summary is 
given. 



II. CONTINUOUS TIME RANDOM WALK 

In the decoupled version of the continuous time ran- 
dom walk (CTRW), a random walker hops from site to 
site and at each site it is trapped for a random time . 
For this well known model two independent probability 
densities describe the random walk. The first is 4'{'t)i the 
probability density function (PDF) of the independent 
identically distributed (IID) pausing times between suc- 
cessive steps. The second is the PDF / (r) for the IID 
displacements of the random walker at each step. Thus 
the CTRW describes a process for which the particle is 
trapped on the origin for time ti, it then jumps to ri, it 
is trapped on ri for time T2 and then it jumps to a new 
location, the process is then renewed. In what follows 
we assume /(r) has a finite variance and a zero mean. 
The asymptotic behavior of the decoupled CTRW is well 
investigated P,p|, |3l| , p4| -p^ and we now summarize some 
results from the CTRW literature which are of relevance 
to our work. 

Let P (r, t) be the PDF of finding the CTRW particle 
at r at time t. Let Nct {s, t) be the probability that s 
steps are made in the time interval (0, t) so clearly 



Y,NcT{s,t) = l 



(3) 



s=0 



and the subscript CT denoted the CTRW. Because the 
model is decoupled 



P{r,t)^Y.^CT{s,t)W{r, s), 



(4) 



and W (r, s) is the probability density that a particle has 
reached r after s steps. In what follows we shall use the 
Fourier (r — > k) and Laplace {t u) transforms, we use 
the convention that the argument of a function indicates 
in which space the function is defined, e.g.. 



P{r,u)^J2NcT{s,u)W{r,s) 



(5) 



s=0 



is the Laplace transform of P(r, t). 

W (r, s) will generally depend on /(r), however we are 
interested only in the large time behavior of P (r, t) mean- 
ing that only contributions from large s are important. 
From standard central limit theorem we know 

I^(r,.)^,^^G(r,.) = ^^-i^exp(^-y . (6) 

We have used the assumption that the system is unbi- 
ased and use convenient units. For most systems and for 
large times NcT{s,t) is concentrated on s = t/rav and 
Tav — Jq t'4'{t)dt is the averaged pausing time. In this 
case P{r,t) will become Gaussian when t ^ 00. When 
Nct (s, t) is broad a non-Gaussian behavior is found. 
This is the case when Tav = 00 or in other words when 
the Laplace transform of behaves as 



V'(m) - 1 -m" -f •• • < a < 1, 



(7) 



and we have used convenient units. Shlesinger showed 
that in this case an anomalous sub-diffusive behavior is 
found, (r^) ~ t". Because the steps are independent and 
based on convolution theorem of Laplace transform 



Nct {s,u) 



1 - 1/' (w) 



exp {s In [■)/) (u)W 



exp (— su") , 



(8) 



and Nct{s,u) is the Laplace transform of NcT{s,t)- 

Following we replace the summation in Eq. (||) with 
integration and use Eq. m) to find 



(9) 



P{r,u) ^ I n{s,u)G{r, s)ds 
Jo 

and according to Eq. (||) 

n{s, u) = u"'^ exp (— su") . 
Using the inverse Laplace transform we find 

/•oo 

P{r,t) - / n{s,t)G{r,s)ds 



(10) 



(11) 



and n{s, t) is the inverse Laplace transform of n(s, u) 

iHil 



n(s, t) 



1 t 



a s 



1+1/c 



5I/C 



(12) 



la {z) in Eq. ( |12[ ) is a one sided Levy stable probability 
density whose Laplace transform is 



la (u) ^ e ""^/q {x) dx ^ e 
Jo 



(13) 
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According to Eq. (|Tl|) the large time behavior of the 
CTRW solution is reached using an integral transforma- 
tion of the Gaussian solution of ordinary diffusion pro- 
cesses [i.e., of G(r, s)]. As we shall see similar transfor- 
mations can be used to solve the FFPE and in particular 
the fractional diffusion equation. 

The kernel n (s, t) is a non negative PDF normalized 
according to 



n (s, t) ds = 1, 



(14) 



it replaces the CTRW probability Nct{s, t). Notice that 
the PDF n{s,t), like Nct{s, t), is independent of the di- 
mensionality of the problem d. Some properties of n(s, t) 
are given in All moments of n(s,t) are finite and 

are given later in Eq. (|^. It is easy to show that 
n{s,t) = [dN{s,t)/ds] with N{s,t) defined in Eq. (|). 
n(s, t) is called the "inverse" one sided Levy stable prob- 
ability density | |37[ |. 

When a = 1, we find n{s, t) = 6{s—t) and hence P(r, t) 
is Gaussian. This is expected since when a = 1 the first 
moment of pausing times Tav is finite, therefore the law 
of large numbers is valid, and hence we expect that the 
number of steps s in the random walk scheme will follow 
s ^ t/Tav When a < 1 the law of large numbers is 
not valid and instead the random number of steps s is 



described by ri(s, t). Thus the transformation, Eq. (11), 
has a meaning of a generalized law of large numbers. 

Eq. ( |l^ ) gives the kernel n{s, t) in terms of a one sided 
Levy stable density. Schneider has expressed Levy 
stable densities in terms of a Fox H function |4ll.[42| 



la{z) 



1 Tjlfi 



(-1,1) 
(-l/a,l/a) 



(15) 



Asymptotic behaviors of the one sided Levy density can 
be found in Feller's book or based on the asymp- 
totic behaviors of the H Fox function ||^,|l|,|2|. For the 
cases a = 1/3, 1/2, 2/3, closed form equations in terms of 
known functions, may be found in |39| (notice that ||36| ] 
points out relevant errors in the literature on Levy stable 
densities). 

Important for our purposes is the result obtained by 
Tunaley j34| already in 1974 which expressed the asymp- 
totic behavior of the CTRW solution P (r, t), in terms of 
its Fourier-Laplace transform, as shown in ]34| 



P(k,«) 



(16) 



In Appendix A we verify that Eq. (|T^) is indeed the 
Fourier-Laplace transform of Eq. (|l^) . The inversion of 
equation (|l^) was accomplished by Tunaley in one 
dimension and by Schneider and Wyss p6|] in dimensions 
two and three (see more details below). 



A. CTRW Solution 

Let us consider an example of the (decoupled) CTRW 
process. The solution of the CTRW in k, u space is M 



1 - iIj{u) 



1 



(17) 



Usually CTRW solutions are found based on numerical 
inverse Fourier-Laplace transform of Eq. ( [l7| ) . Here we 
find an exact solution of the CTRW process for a spe- 
cial choice of f{r) and '0(^)- Our solution is an infinite 
sum of well known functions. We assume the PDF of 
jump times ip{t) to be one sided Levy stable density with 
'0(u) = exp(— u") Displacements are assumed to be 
Gaussian and then 



W (r, s) 



{ins 



(18) 



is exact, not only asymptotic. For this choice of PDFs 
the solution of the CTRW can be found explicitly. We 
use 



NcTis,u) 



1 — cxp(— u") 



exp {—su") 



(19) 



and the convolution theorem of Laplace transform to find 



Nct{s, t) = Lc 



Noit) 
t 



1 



,1/a 



{t) 

— La 



(s + l)l/c 



and 



La. (t) 



la{t)dt 



(20) 



(21) 



is the one sided Levy stable distribution. Inserting Eqs. 
(|l|), (|0l) in Eq. (|) we find 



P{T,t) = [1-La {t)]S{r) + 



51/Q 



La 



(s+1) 



1/q 



{Atts 



(22) 



The first term on the right hand side describes random 
walks for which the particle did not leave the origin 
within the observation time t, the other terms describe 
random walks where the number of steps is s. In Fig. |l| 
we show the solution of the CTRW process in a scaling 
form. We consider a three dimensional case, a = 1/2 and 
use 
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Li/2 (i) = ^r[l/2,l/(4t)] 



(23) 



Here r(-,-) denotes the incomplete Gamma function. 
The figure shows r^P{r,t) versus the scaling variable 
f — r'^ /t". For large times t the solution converges to the 
asymptotic solution found based on the integral transfor- 
mation Eq. (|rT|). 




FIG. 1. We show r^P{r,t) versus C = /t" for the CTRW 
process and for a = 1/2, d = 3. The curves in the figure, t — 2 
(dot dashed), t = 20 (dotted), t = 200 (dashed) and t = 2000 
(stars) are converging in the limit t ^ cx) to the asymptotic 
master curve (solid) . This master curve was computed numer- 
ically based on Eq. ( pd| ) . It is difficult to distinguish between 
the results for t = 200, 2000 and the asymptotic result. Not 
shown is the delta function contribution at r = 0. 



Crn.dTil + m) 



T{1 + am) ' 



(27) 



To derive Eq. ( p7| ) we used the small u expansion of the 
Laplace transform of Eq. ( p5|) and Tauberian theorem. 
Later we shall show that the moments in Eq. ( p7|) are 
identical to the moments obtained directly from the inte- 
gral transformation Eq. (pH). Hence our interpretation 
of Eq. ( |Tl| ) as the asymptotic solution of the CTRW is 
justified, [however our derivation of Eq. (|2^), is based on 
a specific choice of 2p{t) and /(r)]. 



III. FRACTIONAL DIFFUSION EQUATION 

The fractional diffusion equation describes the asymp- 
totic behavior of the CTRW Eq. (|l^). Balakrishnan 
p3| has derived a fractional diffusion process based upon 
a generalization of Brownian motion in one dimension. 
Schneider and Wyss ||2^,0 have formulated the follow- 
ing fractional diffusion equation describing this process 



dP (r, t) 
dt 



=0 Dl-'^W'P{r,t), 



(28) 



and < a < 1. The fractional Riemann Lioville deriva- 
tive oDl-°' in Eq. is defined [|| 



oDl-"Z{t) 



1 8 r^^, Z(t') 



T{a)dtJo {t-t'y- 



(29) 



The Fourier-Laplace solution of the fractional diffusion 
equation, with initial conditions P{r, t = 0) = 5{r), 



Let us calculate the Cartesian moments 
M(2mi,---,2mrf) = 

poo 

dxi • • • / dxdxl""' ■ ■ ■ xf'^Pixi, ■■■,xd,t) (24) 



with non negative integers mi,TO2,---. Clearly the 
odd moments are equal zero and from normalization 
Af(0, • • • , 0) = 1. In Appendix B we find 



Cm.d ^ ^ \ Lq 



M (2mi,---,2md) 



+ 



1/a 



with 



Cm,d 



(25) 



(26) 



and m — nii > 0. In the Appendix we also show that 
for t — > oo 

M(2mi,---,2md) ~ 



P(k,u) = 



k2 



(30) 



is identical to the right hand side of Eq. ( p^ ) . Therefore 
the integral solution of fractional diffusion equation is 



P(r,t) 



n(s, t)G (r, s) ds. 



(31) 



Thus solution Eq. (|11|) which was only asymptotic in the 
CTRW framework is exact within the fractional diffusion 
equation approach. The integral solution Eq. ( ^ ) for 
d = 1 was investigated by Saichev and Zaslavsky, we 
investigate also the cases d = 2, 3 which exhibit behaviors 
different than the one dimensional case. 

As mentioned, Schneider and Wyss [p6| have inverted 
Eq. ( ^0| ) finding the solution to the fractional diffusion 
equation in terms of rather formal Fox function. Later 
we shall show that moments of P{r, t) calculated based 
on the integral solution Eq. ( pl| ) in dimensions d— 1,2,3 
are identical to the moments calculated based upon the 
Fox function solution. In this sense we show that the two 
solutions are identical. 

At this stage we have only shown that the two ap- 
proaches, asymptotic CTRW and fractional diffusion 
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equation are identical. The reader might be wondering 
why should we bother with the fractional equation if the 
CTRW approach gives identical results. The situation 
is some what similar to the standard diffusion equation 
which predicts a Gaussian evolution. The diffusion equa- 
tion a = 1 is an asymptotic equation describing much 
more general random walks. The main advantage of the 
diffusion equation over a random walk approach is its 
simplicity. Solutions with special boundary conditions 
(reflecting/absorbing) are relatively simple and usually 
capture the essence of the more complex random walks. 
Another extension of the diffusion equation is the diffu- 
sion in external field as described by Fokker-Planck equa- 
tions, such an extension for random walks is cumbersome. 
The same is true for the fractional diffusion equation. It 
can serve as a phenomenological tool describing anoma- 
lous diffusion. As we shall show in the next section we 
can include the effects of an external field. On the other 
hand CTRW by definition is not built to consider an ex- 
ternal field. 



A. Integral Solution 

We first notice that the integral solution Eq. ( |3l| ) 
shows that P (r, t) is normalized and non-negative. The 
normalization is easily seen with the help of Eq. ( p^ ) 
and the non negativity of P (r, t) is evident because both 
n(s,t) and G{r,s) are non negative. 

Calculation of spherical moments {r"{t)) is now con- 
sidered. The calculation follows two steps, the first is to 
calculate {r^{s)) for Gaussian diffusion (or find the Gaus- 
sian moments in a text book) then using the transforma- 
tion defined in Eq. ( |3l] ) we find the moments {r^{t)) for 
the fractional case. More precisely in Laplace u space 



(r" (u)) = J dsn{s, ") y G (r, s) r"dr 
the Gaussian spherical moments are 



(32) 



G (r, s) r"dr = ildT [ ^^^^ ) 2"s"/2 (33) 



with ill = I/a/tTj = 1 and Q3 — 2/^/Tr. The moments 
for the fractional case are found using Eq. ( ^ 

(r" (u)) = 



J7^2"r ( adL^] r(l + n/2)w-i-""/2. (34) 



Using the inverse Laplace transformation we find 



r(^)2"r(i + n/2) 

r(l + Q!n/2) 



For n = 2 



ir' {t))= 



2dt° 



r(i + a)' 



(35) 



(36) 



When a = 1 the moments in Eq. ( p5| ) are identical to 
the Gaussian moments in Eq. (|3^). 

In Appendix C we calculate the Cartesian moments 
defined in Eq. (|2|). We find 



(37) 



V (\ A-m') 

M (2mi, • • • , 2md) = , V "-" 

1(1-1- am) 



which is identical to the right hand side of equation (|27 
Thus the moments of the integral solution of fractional 
diffusion equation arc identical to the asymptotic behav- 
ior of the moments of the CTRW found in the previous 
section. 

Eq. (|3^) was also found in j2^ based on Fox function 
solution of fractional diffusion equation (see details next 
subsection). Since the integral solution Eq. ( ^l| ) gives 
identical moments to those found using the Fox function 
solution one can assume that the two solutions are identi- 
cal. It would be interesting to show in a more direct way 
that the Fox function solution is identical to the integral 
solution. In Appendix D we use a theorem on Fox func- 
tions and attempt to give such a proof. Unfortunately 
we fail. The interested reader is referred to Appendix D. 



B. Fox Function Solution 

In [^, the Mellin transform was used to find the so- 
lution of the fractional diffusion equation in terms of Fox 
H function 

P{r,t) =a-^TT-'^/'^r-'^x 

I (d/2,l/a),(l,l/a); 
The asymptotic expression for this solution is |^ 

P (r, t) ~ K°'r-'^£_m^) exp (^-Ai^^) , 
where f = r^/P is the scaling variable, 

= ^-dn^-^ (2 - a)-'/' ^[«(d+l)/2-l]/(2-a)] (40) 

and 

Ai = (2 ~ a) a"/(2-«)2-2/(2-a). (41) 



(38) 



(39) 



Eq. (39) is valid for ^ >> 1. For a brief introduction to 
Fox functions and its application see Gel. 
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The behavior for ^ << 1 is found in Appendix E based 
on the asymptotic expansion of the H function | ^ , ^ . 
For dimensions d — \ we find 



(_l)»^"/2 



n— /I i 



(42) 



and for d = 3 



P (r, t) 



_^ oo 



47rt3a/2£i/2 „!r [1 - a(l + nl2)\ ' 

^ n— / /J 



(43) 



The leading terms in these expansions are for d — \ 

1 1 



P(r,t) 
and for c? = 3 

P(r,t) = 



2r(l-a/2) <"/2 



47rr(l-a) 



(44) 



(45) 



We see that for d = 1 and a = 1, P(r = 0, t) = l/(2V7rt), 
as expected for this normal case. We also see that for 
d = 3, a 7^ 1 and when r ^ the solution diverges like 
P(r, t) ~ l/r. This behavior is not unphysical and P(r, t) 
is normalized according to 47r J P(r, t)r^dr = 1. 

For d = 2 the asymptotic expansion of the Fox function 
is not known (see more details in Appendix E). Using Eq. 
( pl| ) one can show that for d = 2, and ^ << 1 



P(r,t) 



1 



7rr(l - Q;)t" 



In 



(46) 



) were derived independently by A. I. Saichev 



Eqs. 
11- 

The CTRW behavior on the origin is different than 
what we have found for the fractional diffusion equation 
Eqs. (|§0). Within the CTRW on the decay on the 
origin is described with the first term on the right hand 
side of Eq. (|22|), [I - (t)] S (r). This term describes 
random walks for which the particle is trapped on its 
initial location during the time of observation t. This 
CTRW term decays like t^"6 (r) for long times and no 
such singular term is found in the solution of the frac- 
tional. In dimensions d = 2, 3 the decay of the CTRW 
singular term is as slow as the decay of P(0,t) found 
from the fractional diffusion equation. Hence on the ori- 
gin the fractional diffusion approximation does not work 
well. In contrast, for normal random walks, the singu- 
lar term decays exponentially with time and then the 
diffusion approximation works well already after an ex- 
ponentially short time. 



C. Example 



solution are not tabulated and hence from a practical 
point of view one has to use numerical methods to find 
the solution. Also the convergence of the asymptotic so- 
lution for ^ << 1 is not expected to be fast and hence 
the integral solution is of special practical use. 
We consider the case a = 1/2 and d = 3. Using 



^1/2 (z) 



1 



z^^/^exp [-l/(4z)] 



(47) 



with z > 0. We find the integral solution P(r, t) Eq. ( |3l| ) 
using numerical integration. We have used Mathematica 
which gave all the numerical results without difficulty. In 
Fig. H we show P(r, t) versus r on a semi log plot and for 
different times t. Close to the origin r = we observe a 
sharp increase of P(r, t) as predicted in Eq. (^5|). 

More detailed behavior of P(r, t) is presented in Fig. 
^ where we show rP{r^ t) versus r. In Fig. ^ we also 
exhibit linear curves based on the asymptotic expansion 
Eq. (O) which predicts 



rPir,t) Ci{t) ~ rC2{t) 



where 



Ci{t) = 



1 



(48) 



(49) 



and C2(<) = l/[47rr(l/4)t3/4]. Eq. (||) is valid when 
^ = r^/t^^^ « 1 and as expected we see in Fig. ^ 
that for this case the numerical integral solution and the 
asymptotic solution (Esh agree well. 



FIG. 2. P{r,t) versus r on a semi log plot and for different 
times t = 0.02 (dotted), t = 2 (dash dotted) and t = 200 
(dashed). Note that P{r — 0,t) — oo hence the solution in 
figure is cutofi^ close to r — > 0. 



Here we use the integral solution of the fractional dif- 
fusion equation for a specific example. The Fox function 
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FIG. 4. Scaling behavior of P{r,t) for a = 1/2 and d = 3. 
We show r^P{r,t) versus ^ = r^/t" for times t = 0.02 (dot- 
ted), t — 2 (dash dotted) and t = 200 (dashed), due to exact 
scaling of the solution one cannot distinguish between the 
three curves in the figure. Also shown (solid curves) are the 
two asymptotic behaviors valid for small and large values of 
the scaling variable ^. For ^ >> 1 we used Eq. ( p9| ) and for 
^ << 1 we used the leading term in the expansion Eq. (^3|). 




FIG. 3. Behavior of P{r,t) close to the origin r = for 
a = 1/2 and d = 3. We show rP{r,t) versus r for different 
times t = 0.02,0.2,2,200 (dashed curves). Because P{r,t) 
diverges when r ^ like P{r,t) ~ 1/r, the figures shows 
that rP{r,t) — > C\ (i). Also shown (linear curves) the ap- 
proximation Eq. (ksl), as expected when ^ = << 1 the 
approximation is in good agreement with the integral solu- 
tion. 

Finally in Fig. ^ we show our results in a scaling form, 
similar to the way we presented the CTRW results in 
Fig. which shows the CTRW solution for both short 
and long times. We present r^P{r, t) versus ^ for different 
choices of time t. We observe collapse of all curves, both 
for shorter and longer times, onto one master curve. We 
also show the asymptotic behaviors ^ >> 1 and ^ << 
1, Eq. ( |39| ) and Eq. (Esj) respectively. The numerical 
integral solution Eq. (|3lD agrees well with the asymptotic 
behaviors in the appropriate regimes. Comparing Fig. 
(|l|) and Fig. (||) we see that the solution of the fractional 
diffusion equation approximates the exact solution of the 
CTRW very well when i — oo and r ^ 0. 
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IV. FRACTIONAL FOKKER-PLANCK 
EQUATION 

Let us consider the fractional Fokker-Planck equation 
(FFPE) |3^] describing the stochastic evolution of a test 
particle under combined influence of external force field 
F{x) and a thermal heat bath. The equation reads 



and 



d F{x) d'^ 
kbT ~^ 



(50) 



(51) 



is the Fokker-Planck operator. Ka is a generalized dif- 
fusion constant and T is temperature. The fractional 
derivative qD]'" was defined in Eq. (p9|). We consider 
the one dimensional case and extensions to higher dimen- 
sions are straight forward. When F{x) — the FFPE 
Eq. ( pO| ) reduces to the fractional diffusion equation (|2^). 
When a — 1 we recover the ordinary Fokker-Planck equa- 
tion. The stationary solution of the FFPE is the Boltz- 
mann distribution and the equation is compatible with 
linear response theory | ]30| |. 

In this section we find an integral solution of the FFPE. 
We show that the solution is normalized and non neg- 
ative, an issue not discussed in Without loss of 
generality we set Ka — 1. We consider initial condition 
P{x,t = 0) =S{x-xo). 

First let us show that the solution is normalized. Inte- 
grating Eq. ( [50| ) with respect to x and using the bound- 
ary conditions 



P(x,i)U=±oo = ^^^U = ±oo = 



we find 



P{x,t)dx = 0, 



(52) 



(53) 



meaning that normalization is conserved as it should. 

We write the solution of the FFPE Eq. (|5^) in terms 
of an integral of a product of two functions 



/•oo 

P{x,t) = i dsfi{s,t)Pi 
Jo 



{x,s) 



(54) 



where 
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dPi{x,s) 
ds 



= LppPi (x, s) . 



(55) 



Piix, s) is a normalized solution of the ordinary Fokker- 
Planck equation with initial conditions Pi{x,s = 0) = 
6{x — xq). Methods of solution of Eq. ( |55| ) are given in 
Risken's book |47|. In what follows we shall prove that 
n{s,t) = n{s,t), defined in Eq. jl^). 

We use the Laplace transform of Eq. (^4|) and normal- 
ization condition of P{x, t) to show 



7T,(s, U)ds — \/u. 



(56) 



Hence ri(s, t) is normalized according to n{s, t)ds = 1. 
The Laplace transform of Eq. (M) reads 



(57) 



uP{x,u) — 5{x — Xq) — °'LfpP{x,U), 

inserting Eq. ( |5^ ) in Eq. (|57| ) we find 



u I n{s,u) Pi {x,s) ds — S{x — xq) 
lo 



n{s, u)LfpPi (x, s) ds, 
integrating by parts using Eq. (|5^) we find 

n (s, u) Pi {x, s) ds — S{x — xq) 



(58) 



" [h (oo, u) Pi {x, s — oo) — h (0, u) Pi {x, 0)] 



,1-a 



—n{s,u) 

OS 



Pi(x, s)ds. 



(59) 



From Eq. ([Sq ) h{oo, u) = and since Pi{x, 0) = S{x—xo) 
we may rewrite Eq. (|59|) 



Uh{s, u) + M"'" 



— n (s,-u) 
as 



Pi(x, s)ds 



[l-Mi^"n(0,u)] (5(x-a;o). 



(60) 



Eq. ( |60D is solved once both of its sides are equal zero; 
therefore two conditions must be satisfied, the first 



h{0,u) = 



and the second 



d 

—h (s, u) — —u"n (s, u) . 



(61) 



(62) 



The solution of Eq. ( |62| ) with initial condition Eq. ( |6l| ) 
is 



n(s,M) = u" "'^ exp (— w"s) . 



(63) 



Thus n(s, u) = n(s, u) found in the context of the solution 
of CTRW Eq. ®. 

We see that the integral solution of the FFPE has a 
similar structure as the solutions of the fractional diffu- 
sion equation or equivalently of the asymptotic CTRW. 
The solution is Eq. (||) with n(s, t) defined in Eq. (|l|) 
and Pi{x, s) being the solution of corresponding ordinary 
Fokker-Planck equation. It is easy to understand that the 
solution, P{x, t) is normalized and non negative because 
both Pi{x,t) and n{s,t) are normalized PDFs. 

In a different method of solution of the fractional 
Fokker-Planck equation was investigated. It can be 
shown that P{x, t) can be expanded in an eigen function 
expansion, which is similar to the standard eigen function 
expansion of solution of ordinary Fokker-Planck equa- 
tion, however eigen modes decay according to Mittag- 
LefBer relaxation instead of the standard exponential de- 
cay of the modes in the integer Fokker-Planck equation. 
It can be shown that the integral solution investigated 
here is identical to the eigen function expansion in |30[| . 



A. Example 1, Biased Fractional Wiener Process 

Consider the biased fractional diffusion process defined 
with a generalized diffusion coefficient Ka and a uniform 
force field F{x) = F > 0. For this case the mean dis- 
placement grows slower than linear with time according 
to 



{x{t)) 



(64) 



kbTr{l + a) 

The well known solution of the ordinary Fokker-Planck 



equation is 



Pi (x, s) = - — exp 
47rs 



X - Fs 



As 



(65) 



Eq. ( |65| ) describes a biased Wiener process and F = 
F/{k\jT). The solution P{x,t) for the fractional case is 
found using the transformation Eq. (|5^). In Laplace u 
space the solution is 



P{x,u) 



^1 + 4(ut)' 



: exp 



F{x - ^1 + 4(ur)"|a;|) 



(66) 



and r" = l/{FKa)'^. For (rit)" << 1, corresponding to 
the long time behavior of the solution P{x,t), we find 



P{x,u) 



Fu°'-'^t" exp \^-Ft"u"x 
Fu°'-^T°' exp (-F|a;| - F|a;|r"M" 



x>0 

x<0. 
(67) 
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Since P{x, u) ^ u°' ^ and hence according to Taube- 

rian theorem P{x,t) ^ l/t"' we have P{x,t) = for 
a; < when t — > oo. Therefore, using the inverse Laplace 
transform of Eq. (|67[) the asymptotic behavior of P{x, t) 



Hm P{x,t) 



1 t 1 ( t \ 



< a; 
> 



(68) 



and j4 = _Fr" . Integrating Eq. ( pq ) we find the distribu- 
tion function 



hm 

t — >oo 



P(a;, i)(ia; = 1 — 



(69) 



valid for a; > 0. Eq. ( p^ ) was derived also in ||3^ based 
on the biased CTRW, thus as expected the solution of 
the fractional Fokker-Planck equation converges to the 
solution of the CTRW in the limit of large t. 

On the origin one can use Tauberian theorem to show 



P(0,0 



A 



ni-a) 



t~ 



(70) 



valid for long times. For the case = we have found 
P(0, t) ~ so as expected the decay on the origin is 

faster for the biased case since particles are drifting away 
from the origin. 

In Fig. ||we present the solution for the case a = 1/2, 
then 



P(x,t) 



ds- 



: exp 



{x ~ Fsf 



4s 



(71) 



which is evaluated numerically. For large times we have 



P(x,t) 



A 



exp 



A^ 



(72) 



for < x. As seen in the figure the exact result exhibits a 
strong sensitivity on initial condition and the maximum 
of P(a;, t) is located on a; = 0. This is different than ordi- 
nary diffusion process in which the maximum of P(a^, t) 
is on (a;(t)). The curves in Fig 
observed by Scher and MontroU 
ulation of CTRW and also by Weissman et al who 
investigated biased CTRW using an analytical approach. 
The FFPE solution presented here is much simpler than 
the CTRW solution, still it captures all the important 
features of the more complex CTRW result. 



are similar to those 
based on lattice sim- 
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FIG. 5. The dynamics of P{x, t) for the fractional diffu- 
sion process in external field F with a — 1/2. We show 
the solution for four different times. The dotted curve is 
the asymptotic solution valid for large t. The maximum of 
P[x,t) is for all times on the initial condition a; = 0. We use 



F ^kbT = Ki 



/2 



1. 



B. Example 2, the Fractional Ornstein-Uhlenbeck 
Process 

We consider as a second example the fractional 
Ornstein-Uhlenbeck (OU) process, namely the motion 
of the test particle in harmonic oscillator. This case can- 
not be analyzed using the CTRW in a direct way. The 
CTRW formalism considers only uniformly biased ran- 
dom walks and the fractional processes in non-uniform 
fields are not uniformly biased. We consider a = 1/2, 
F{x)/khT — —X and use the well known solution of the 
ordinary OU process 



Pi [x, s) 



1 



-2s 



exp 



{x-xoe 



2(1 



-2s 



(73) 



The solution of the fractional OU process is then found 
using numerical integration of Eq. (|5^) using Eqs. ( |47| ) 
and ([73|). Our results are presented in Fig. ^. We have 
considered an initial condition xq — 1/2 and we observe 
a strong dependence of the solution on the initial con- 
dition. A cusp on a; = a;o is observed for all times t, 
thus the initial condition has a strong influence on the 
solution. The solution approaches the stationary Gaus- 
sian shape slowly in a power law way and the solution 
deviates from Gaussian for any finite time. Unlike the 
ordinary Gaussian OU process, the maximum of P(a;, t) 
is not on the average {x{t)) but rather the maximum is 
for short times located on the initial condition. 
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FIG. 6. The dynamics of P{x,t) for the fractional Orn- 
stein-Uhlenbeck process and ct = 1/2. We show the solution 
for four different times. The dotted curve is the stationary 
solution, i.e. the Boltzmann distribution which for the Har- 
monic potential is Gaussian. Notice the cusp on the initial 
condition x = xq = 1/2. 

Like the ordinary OU process the fractional OU pro- 
cess has a special role. The ordinary OU process de- 
scribes two types of behaviors, the first is an over damped 
motion of a particle in harmonic potential, the second 
is the velocity of a Brownian particle modeled by the 
Langevin equation, the latter is the basis for the Kramers 
equation. In a similar way the fractional OU describes an 
over damped and anomalous motion in Harmonic poten- 
tial considered in this section and in it can also be 
used to model the velocity of a particle exhibiting a Levy 
walk type of motion |Q. The fractional OU process is 
the basis of the fractional Kramers equation introduced 
recently by Barkai and Silbey |Q , this equation describes 
super diffusion while in this work we have considered sub 
diffusion. 

V. SUMMARY 



and it can serve as a practical tool for finding solution of 
certain fractional kinetic equations. 
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VII. APPENDIX A 

We rewrite Eq. ( ]T^ ) 

/•oo 

P(k,u)-M"-M e-n"°+k')ds, (74) 

^0 

the inverse Fourier transform of Eq. ( |7^ ) is 

P (r, u) - u"-^T-^ |^°° e-"("°+''')ds| , (75) 

and !F~^ is the inverse Fourier transform. Changing the 
order of integration over parameter s and the operation 
(this is later justified by the identity of moments of 
the integral solution and that found in ||26| ) we find using 
Eq. (I) 

P(r,u)^ / n(s,u) — ■ exp (— r^/4s) ds. (76) 

Jo {AttsP^ 

Applying the inverse Laplace transform to this equation, 
changing the order of integrations over s and the inverse 
Laplace operation, we find Eq. (p]). 

VIII. APPENDIX B 

We investigate the Cartesian moments of the CTRW. 
These moments are 

M(2mi,---,2md) = 



Fractional diffusion equation is an asymptotic equa- 
tion which predicts the behavior of the decoupled con- 
tinuous time random walk in the sub diffusive regime. 
The fractional Fokker-Planck equation considers such a 
sub-diffusive motion in an external force field and close 
to thermal equilibrium. In this work we have considered 
an integral transformation which gives the solution of the 
fractional Fokker-Planck equation in terms of solution of 
ordinary Fokker-Planck equation. Solution of ordinary 
Fokker-Planck equation can be found based on different 
analytical and numerical methods p^-|5l[|. The integral 
transformation describes also the long time behavior of 
the CTRW in dimension d = 1,2,3. Thus the transfor- 
mation maps Gaussian diffusion onto fractional diffusion, 



dxi 



2mi m2 ^2m 



51/Q 



t 



(47rs) 



d/2 



cxp 



+ 

2 



1/q 



Xt + ■ ■ ■ X' 



As 



(77) 



Changing order of integration and summation, we use 
the identity 



dxi ■ 



dxr, 



1 



' (47rs) 



d/2 



exp 



4s 
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. . . ^2md _ (-1 r. 



were Cm.d is defined in Eq. (^6|) Inserting Eq. (| 
Eq. we find Eq. @. 

The Laplace transform of Eq. is 



(78) 
^ in 



1 -e- 



M(2mi,---,2md) = C,n,d 
We use 

oo / , s m oo 

^e-'*"^™ = (-1)™ ( — ) XI '=""''1^ 



Xe-'^^^s". (79) 



(-1)"' iz 



da; / 1 — e ^ 

and for small we find 

M (2toi, • • • , 2TO<i) ~ C„,dr(m + l)u-i-"". 

Applying Tauberian theorem [i.e., inverting Eq. 
findEq. (|^. 

IX. APPENDIX C 

We consider the moments 

M (2mi,---,2m<i) = 



(80) 



(P)] we 



dXrl 



n (s, i) G (r, s) 



. . . 2;^'" 



(82) 



changing the order of integration over s and the d dimen- 
sional integration, using Eq. ([78|) we find 



M (2mi,---,2m<j) = 
C„i4 / s^n{s,t)ds. 



(83) 



were Cttj ,d is defined in Eq. (|2^) and m = "ij, 
Using Eq. (|l|) the integral in Eq. (H) is 



OO 



J™(t) = / s'^n{s,t)ds 



1 / t 



I \S I 







(84) 



Notice that /m(0 i^ the m moment of n{s,t). Changing 
the integration variable, in Eq. (84), according to y = 
i/s^/" we find 



^0 



(85) 



the calculation of the negative moment in Eq. ( ^ ) is 
straight forward using Laplace transform technique [^j, 
we find 



Irn{t) 



j-am poo 



T (am) Jq 



u^'^-Ha (u) du 



(86) 



with la (u) = exp(— u") being the Laplace transform of 
the one sided Levy density. Therefor we find 



rji + m) 

r (1 + am) 



(87) 



Imit) 



Inserting Eq. (^Tp in Eq. (g3|) we find the result Eq. (|37 
Notice that when a — 1, Im{t) = t™, this is expected 
since for this case n(s, t) = 5{s — t). 



X. APPENDIX D 

Integral formulas involving the product of two H Fox 
functions are a helpful tool with which H functions can 
be represented in terms of known functions. According 
to El equation (5.1.2) 



/ dx{x^~'H;^^- 

Jo 



xH 



(cj'7i)i^ 
{dj,5j)^ 



] = 



-nTTm+M,n+N r a i 
* ^p+P,q+Q I I 

provided that seven conditions are satisfied, cr > 0, 



Q P 



(88) 



(89) 



(90) 



^ = E"j - E "j+E/^^- E ^j >0' (91) 

J = l j=n+l j=l j=m+l 



M Q N P 

^'-E'^^- E '5.+E^^- E ^^>0' (92) 

i=l j = A/+l 'j = l j=N+l 
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max [{ttj — 1) /aj] 



l<j<n 



min [dj/Sj] < rj 



(93) 



and 



r] < a min [bj/Pj] ~ max [(cj — 1) /7j] . 
l<j<m l<j<iV 

(94) 

We have considered the case when all parameters are real, 
for the case when parameters may become complex see 

We express the integral solution Eq. (|3l|), in terms of 
an integral of a product of two Fox H functions, after 
some rearrangements and with the use of Eq. (|l^) and 
the following identity 



(95) 



we find 



P(r,i) 



/ n{s,t)G {r, s) ds = 
Jo 



-1 -d/2 -d 



/ dxx-^-^'^Hli[x~^' 
Jo 



(-1,1) 

(-l/a,l/a) 



rrl.O 
^0,1 



' (d/2,1) 



(96) 



Using Eq. ( p8|) one can hope to prove that the integral 
solution Eq. (|3l|) and the Fox Function solution Eq. (pF 



are identical, provided of course that all conditions Eqs. 
( p9[j9^) are satisfied. Unfortunately the integral identity 



Eq. (|8J) cannot be used for this aim because condition 
Eq. (94) is not satisfied, inserting the fractional param- 
eters in Eq. we find 



— l/a < — 1/a 



(97) 



which shows that the conditions are not fulfilled. 

It is worthwhile mentioning that all the other condi- 
tions (|89|-p3|) are fulfilled and that it is possible to prove 
that the integral solution and the Fox function solution 
are identical if we replace < with < in Eq. (|9j) . Using 
Eq. (H), the identity 



1,2 



.I/O I (0,1) 

^ I (0, l/a), (d/2, l/a) 



0-2,0 
-"1,2 



a/a I (1,1) 

^ I (d/2, l/a), (1, l/a) 



(98) 



which can be proven based upon the definition of the Fox 
function, and ||l|,|| 



o- TTm,n 



Tj7n,n 



^1 



, (ai,aj)i.p 
' (&.,/3.)i,, 

(aj,o-"i)i,p 

(^J,^/?j)l,5 



(99) 



we find 



P(r,t) = / n{s,t)Gir,s)ds 
Jo 



rr20 / o-2/a 2/q,-1| (1, 1) 

I (d/2, l/a), (1, l/a) 



(100) 



which is the Fox function solution of Eq. (jSg). We em- 
phasize that this equation was not proven in this Ap- 
pendix because condition in Eq. (M) was not satisfied. 



XI. APPENDIX E 



The Fox function is represented as 



rrm,n 

p,<? 



(101) 



(oi, ai) • • • (op, ap) 
(6i,/3i)...(6„/3g) 

and for our choice of parameters 

m — 2,n — 0,p — l,q~2 



ai = 1, ai = 1 

h = d/2, [3i = l/a, 62 = 1, 132^ l/a (102) 

The asymptotic expansion of the H Fox function, for 
< a; << 1, is defined when two conditions are satis- 



fied p 42|. The first is 



(103) 



and for the case 5 = see [ 4l]j4^ . For our case, defined 
by the parameters in Eq. (102), 5 — 2/a~l>{) when 
< a < 1. The second condition is 



(3h {bj + A) ^ I3j {bh + k) (104) 



for 



j^h] j = h = l,---m; A,fc = 0, 1,2,---. (105) 



Using Eq. (102) condition Eq. (105) reads 
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- (1 + A) ^ - {d/2 + k) 
a a 



(106) 



therefore the condition is satisfied for dimensions d ~ 1 
and d = 'd but not for li = 2. When conditions are 
satisfied 



(ai, ai) • • • {up, ap) 
(6i,/3i)---(5,,/3,) 



m oo 



n;'Li,.//.r (bj - Pj^hm) (i - a, + a,a,fe) 

h=l k=0 



where 



(-if [n]^,„+ir(i-6, + /3,a.fe) 



S.h,k = (bh + k) /Ph. 



(107) 



(108) 



Using Eq. (|107|) we find 



H 



1,2 



(1,1) 

(rf/2,l/a),(l,l/a) 



k=0 



r(l (-l)''a;' 
r (1 - Qd/2 - a/c) fc! 



g. r(ci/2-i-fc)(-i)^x"(i+^) ^^pg^ 



/£ = 



r(l-a-Qfc)fc! 



Using some simple manipulations we find our results Eqs. 
(H) and (H). 
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